!=======================================================
! HYDROTHERMAL ALTERATION OF THERMAL DIFFUSIVITY
! Modified from Luc's code, G. Ito 8/7/06
!=======================================================
!-------------------------------------------------------------
subroutine ReadHydro()
!-------------------------------------------------------------
include 'precision.inc'
include 'params.inc'
common /hydroth/ xmaxdepth,xmaxt,xmaxstr,xenhc1,xenhc2,   &
                   xhyd1,xhyd2,xhyd3,xhyd4,xcond1,xcond2,xcond3

open( 9, file='hydrother.inp',status='old',err=2001 )
if_hydro = 1
call AdvanceToNextInputLine( 9 )
read (9,*) xmaxdepth,xmaxt,xmaxstr,xenhc1,xenhc2
call AdvanceToNextInputLine( 9 )
read (9,*) xhyd1,xhyd2,xhyd3,xhyd4
call AdvanceToNextInputLine( 9 )
read (9,*) xcond1,xcond2,xcond3

write(*,*) '------------------ Hydrothermal Diffusivity ---------------------------'
if (xenhc2.lt.xenhc1) then
  write(*,*) 'WARNING:  xenhc2 corrected to be >/= xenhc1'
endif
write(*,*) ' xmaxdepth,   xmaxt,  xmaxstr,   xenhc1,   xenhc2'
write(*,'(5f10.2)') xmaxdepth,xmaxt,xmaxstr,xenhc1,xenhc2
write(*,*) '    xhyd1,    xhyd2,   xhyd3,   xhyd4'
write(*,'(4f10.2)') xhyd1,    xhyd2,   xhyd3,     xhyd4
write(*,*) '   xcond1,   xcond2,    xcond3'
write(*,'(3f10.2)') xcond1,     xcond2,    xcond3
write(*,*) '-----------------------------------------------------------------------'
close( 9 )


return

2001 if_hydro = 0
write(*,*) '>>NO Hydrothermal effects on Diffusivity<<<'
return

end

!-------------------------------------------------------------
function HydroCond(i,j)
!-------------------------------------------------------------
  include 'precision.inc'
  include 'params.inc'
  include 'arrays.inc'
  common /hydroth/ xmaxdepth,xmaxt,xmaxstr,xenhc1,xenhc2, &
                   xhyd1,xhyd2,xhyd3,xhyd4,xcond1,xcond2,xcond3
  pi = 3.14159265358979323846
  
  iph = iphase(i,j,phasez(j,i))
  tmpr = 0.25*(temp(j,i)+temp(j+1,i)+temp(j,i+1)+temp(j+1,i+1))
  yc = 0.25*(cord(j,i,2)+cord(j+1,i,2)+cord(j,i+1,2)+cord(j+1,i+1,2))
  xc = 0.25*(cord(j,i,1)+cord(j+1,i,1)+cord(j,i+1,1)+cord(j+1,i+1,1))
  xc = dabs(xc)
    
  if( tmpr.le.xmaxt .and. yc.ge.xmaxdepth) then
      if (xcond1.lt.0.0) then			!changed 5/07
        cdum=dmin1(aps(j,i)/xmaxstr, 1.0)
        HydroCond=(xenhc1 + cdum*(xenhc2-xenhc1))*conduct(iph) 
      else
        if (xc.le.xhyd1) then
          fdum=xcond1
        elseif (xc.gt.xhyd1.and.xc.lt.xhyd2) then
	  xx=(xc-xhyd1)/(xhyd2-xhyd1)
          fdum=xcond1+(xcond2-xcond1)*(1-dcos(pi*xx))/2.0
	elseif (xc.ge.xhyd2.and.xc.le.xhyd3) then	  
          fdum=xcond2
	elseif (xc.gt.xhyd3.and.xc.lt.xhyd4) then
	  xx=(xc-xhyd3)/(xhyd4-xhyd3)
          fdum=xcond2+(xcond3-xcond2)*(1-dcos(pi*xx))/2.0
	else
	  fdum=xcond3
	endif
	HydroCond=fdum*conduct(iph)
      endif        
  else
      HydroCond=conduct(iph)
  endif
  
!  if (j.eq.1) then
!    write(*,'(2i5,4f15.3)') j, i, xc, xcond1, fdum, HydroCond
!  endif
return

end

